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We study the ground state of the attractive one-dimensional Bose-Hubbard model, and in par- 
ticular the nature of the crossover between the weak interaction and strong interaction regimes for 
finite system sizes. Indicator properties like the gap between the ground and first excited energy 
levels, and the incremental ground-state wavefunction overlaps are used to locate different regimes. 
Using mean-field theory we predict that there are two distinct crossovers connected to spontaneous 
symmetry breaking of the ground state. The first crossover arises in an analysis valid for large 
L with finite A'', where L is the number of lattice sites and A'^ is the total particle number. An 
alternative approach valid for large A'^ with finite L yields a second crossover. For small system 
sizes we numerically investigate the model and observe that there are signatures of both crossovers. 
We compare with exact results from Bethe ansatz methods in several limiting cases to explore the 
validity for these numerical and mean-field schemes. The results indicate that for finite attractive 
systems there are generically three ground-state phases of the model. 



I. INTRODUCTION 

Systems of attractive bosons are one of the most 
intriguing current topics in physics. For instance 
they might lead the way for fabricating mesoscopic 
Schrodinger cat statesii^, and in the experimental con- 
text, they have been used to produce the Bosenova 
phenomenal^. The substantial amount of research un- 
dertaken recently^i^iii^i^i^iii poses new questions sur- 
rounding systems of attractive bosons. An almost ideal 
realisation of a lattice Bose gas - the Bose-Hubbard model 
- has been found in bosons trapped inside optical lat- 
tices'^. The use of techniques like Feshbach resonances 
allows tuning of the scattering length, i.e. changing the 
interaction strengh, even crossing from repulsive to at- 
tractive^. The theoretical boson model predicts a dra- 
matic change in the ground state of a large but finite 
system when the attractive interaction strength is varied 
from weak to strong attractive, see figure [T] and figure [2] 
for visualisation. Currently technical difficulties make ex- 
periments on attractive system considerably harder com- 
pared to repulsive system''^. Once handling and stabil- 
ity of attractive bosons in optical lattices allows exper- 
iments at controlled and varying interaction strengths 
this general transitional feature could be experimental 
verifiable. Potential candidates for measurements might 
include correlation functions, momentum distribution af- 
ter release from the trap and the low-lying energy spec- 
trum. Historically, the theoretical study of attractive 
bosonic systems has received little attention due to dif- 
ficulticaiS fike non-saturation or high site occupancy. A 
number of numerical and approximative studies for a va- 
riety of attractive bosonic systemsii^iiSiiS have found a 



transitional regime between the strong and the weak in- 
teracting regions. This crossover can be seen in prop- 
erties like the energy spectrum^, correlation functions^ 
or entanglement^. All these properties have the common 
feature that the crossover becomes sharper and more pro- 
nounced for larger system sizes, a region where numerical 
and approximative techniques enter a region of uncer- 
tainty. A transition is also seen in studies of mean-field 
techniques of the non-linear Schrodinger equation in the 
context of the Bogolyobov approximation or in solitonic 
solutions of the Gross-Pitaevskii equatio n -'^^i^^ . Despite 
an early Bethe Ansatz solutionis for the continuum Bose 
gas with contact interactions, the exact treatment of at- 
tractive quantum systems lags behind the study of similar 
repulsive systemsHi^iSS. 

In this work we consider the one-dimensional periodic 
Bose-Hubbard model in the attractive regime, as a sim- 
ple boson model with short-range interactions and local 
hopping term. This model is in general not integrable^, 
but it possesses several integrable limits and displays 
rich transitional behavior in the ground state. A quan- 
tum phase transition (QPT) is usually defined as a phase 
transition at T = (i.e. in the ground state) under the 
variation of external parameters, here the attractive in- 
teraction strength. Phase transitions involve taking the 
thermodynamic limit, e.g. having infinitely many par- 
ticles N and lattice sites L. Attractive boson systems 
are conceptually different from repulsive bosons and at- 
tractive/repulsive fermions, in that such a limit cannot 
easily be defined as discussed later on. Nevertheless, for 
large but finite N and L the attractive boson system does 
display an increasingly sharp distinction beween ground- 
state regions, similar to finite size realisations of system. 
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FIG. 1: Generic ground-state behaviour for attractive bosons: 
momentum distribution of trapped bosons, within a mean- 
field approach for L = 50. See also figure [2] for the real space 
"density" and section UlI Al for technical details. In the weak 
interaction region (right front in these figures: r — > 0) the 
system is in an ideal BEG state, all bosons are condensed into 
the lowest momentum and the semi-classical density is flat. 
For strongly interacting bosons (left rear in these figures: r 
1) the momentum distribution is flat, while the translational 
symmetry in real space is broken, all bosons are on the same 
side. In between there is a rich cross over regime which we 
study in this paper. 



We will within this paper denote this generahsation by 
pre-transition and discuss its relevance as a tool for the 
analysis of attractive bosons. 

To characterise the ground-state phases of the model, 
we study two key indicator properties. The first is sim- 
ply the energy gap between the ground and first-excited 
states. For finite systems the gap never vanishes, and 
there is never an occurrence of ground-state broken- 
symmetry in the quantum model as a result. However 
we do observe through numerical analysis that the or- 
der of magnitude of the gap can be significantly different 
across different coupling regimes, which leads to a sense 
of relative quasi-degeneracy^. The second key property 
we study is the incremental ground-state wave function 
overlap, or the fidelity to use the language of quantum in- 
formation theory. Recently there have a been a number 
of papers that have used this concept to study quantum 
phase transitions in the thermodynamic limiij^i^. The 
essence of this approach lies in the fact that if two states 
lie in different quantum phases then they are reliably 
distinguishable, for example through the use of an or- 
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FIG. 2: Real space density for attractive bosons in a trap 
with L = 50 sites in a semi-classical theory. This picture is 
the corresponding density to figure [T] see caption for different 
physical regimes. Technical details text are in section IIII Al 



der parameter. If states are reliably distinguishable then 
they must be orthogonal^ and consequently the wave 
function overlap vanishes. 

For finite systems we propose to modify this approach 
by identifying pre-transitions at couplings for which the 
incremental wavefunction overlap is (locally) minimal, 
see figure 21 For systems which exhibit a quantum phase 
transition in the thermodynamic limit it is then neces- 
sary that the value of the minimum goes to zero in the 
thermodynamic limit. In this manner we can say that 
the occurrence of a minimum in the incremental ground- 
state wavefunction overlap in a finite system is a precur- 
sor for the quantum phase transition in the thermody- 
namic limit. Note that the incremental overlap graphs 
are shown on a unitless axis as the physical interest here 
lies in the existence and location of minima, not the quan- 
titative shape^^. 

The results we find from the study of incremental 
ground-state wavefunction overlaps give overwhelmimg 
support to the mean-field results, viz. the general exis- 
tence of two transitional couplings. Within the context 
of mean-field theory the system exhibits a broken sym- 
metry phase. Our mean-field results point towards the 
existence of two transition couplings, with the critical 
couplings becoming degenerate at zero coupling in the 
limit of large particle number TV and a large number of 
lattice sites L. However by judiciously choosing the scal- 
ing of the parameters our findings also show that the lim- 
its iV ^ oo and L — > oo do not commute. For example 
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the bosonic statistics that underly the system mean that 
it is possible to take iV to infinity while keeping L finite. 
Moreover, we can also take N and L to infinity such that 
the "density" N''-' /L constant for any V > 0. This 
prospect leads to the conclusion that the thermodyamic 
limit of the model appears to not be well-defined. This is 
a significant distinguishing feature compared to fermionic 
lattice systems such as the Hubbard model^^ where the 
thermodynamic limit is well-defined. In section |TT] we in- 
troduce the one-dimensional Bosc-Hubbard model, list 
key properties used for this study and present some nu- 
merical results for small systems. Next a first analysis 
of the pre-transition points is given via different mean- 
field approaches in section HTTl especially the limiting case 
L = 2 and L — > oo are discussed. Section|lV]discusses the 
limiting solvable cases for L = 2 (Bose-Hubbard dimer), 
N = 2 (Haldane-Choy), and L ^ oo (Lieb-Liniger) via 
the Bethe Ansatz solution. The results of the mean-field 
theory and the small size exact diagonalisation are com- 
pared with these exact solutions and the limiting quasi 
root distribution is discussed. The discussion in sectionlVl 
finally puts all three approaches together and concludes 
that in limiting cases, e.g. very small or very large L, 
only two regions might be visible. Nevertheless, our main 
finding is that three ground-state phases exist in the at- 
tractive regime of model for the generic case of finite but 
large number of particles A'^ and lattice sizes L - presum- 
ably the experimentally relevant case. 



pation numbers |ni, 712, n^) with aj\ni, ...rij , ...ul) = 
\Jnj + \ \n\...nj + 1, ...ni). Its dimension d — {N + L — 
1)!/((L — 1)!A^!) grows very rapidly with particle num- 
ber and lattice size. For example, the moderate values 
iV = 10 and L = 20 give the dimension of the Hilbert 
space as d = 20,030,010, strongly limiting exact diag- 
onalisation of systems except for the dimer and trimer 
system. Use of truncation schemes for the dimension and 
quantum Monte Carlo methods arc limited by apriori un- 
known behavior in the transitional regions of interest. 

As the Hamiltonian ([T]) conserves the momen- 
tum^, the matrix representation in a free momen- 
tum basis is block diagonal, and low-lying, quasi- 
degenerate states are characterised by differing mo- 
menta. Defining creation and annihilation opera- 
tors in momentum space via the Fourier transforms 
bk = X]^=i 6xp (2ifcj7r/L)aj, k — 1...L, it can be 

shown these operators satisfy canonical bosonic commu- 
tation relations = Sjkl- The Hamiltonian ([1]), 
acting on a dual lattice of equally L sites (modes) may 
be equivalently expressed as 

^ /27rfc\ 

i7BH= -2i^cos(— j blbk (2) 
fe=i ^ ^ 

L 

— ^ ^ b\b\bmbn5k+l=m+n (modi) 
/c,/,m,n— 1 



II. DEFINITION OF THE MODEL 

We consider a one-dimensional Bose-Hubbard model, 
consisting of bosons with creation (annihilation opera- 
tors) aj (aj) that create (annhiliate) a boson at lattice 
site J, with j running over all L lattice sites. The usual 
bosonic commutation relations such as [Oj, a^] = 5jkl ap- 
ply. For a discussion of the physical origin and limitations 
of this model se e^^'^^ . Particles on the same site interact 
with interaction strength 7. The kinetic term is given by 
nearest-neigbour hopping with coupling strength t, and 
periodic boundary conditions a^+i = ai arc imposed. In 
the real space presentation the Hamiltonian is given by 



The Hamiltonian commutes with the total particle num- 



ber N = X]7=i "-j with Uj 



The physical 



Hilbert space is spanned by Fock states of on-site occu- 



For the remainder of this paper we only consider 
t > and 7 > corresponding to attractive interac- 
tions. The model then incorporates the competition be- 
tween the delocalising and localising effects of the ki- 
netic and the interaction terms respectively In the limit 
7 — the ground state approaches that of non- interacting 
bosons, and is non-degenerate. At the other limit >0 
the ground state becomes L-fold degenerate where the 
ground states consists of N localised bosons on a single 
lattice site, viz. states of the form |0, ...0, rij = N, 0, ...0). 
However for non-zero t the degeneracy is broken and the 
unique ground state is a superposition of these localised 
states, giving rise to a Schrodinger cat state. The low- 
est L energy states in this strong interaction limit form 
a narrow energy band. Within mean-field theory treat- 
ments, as will be shown below, this energy band degener- 
ates at non-zero values of t giving rise to spontaneously 
broken translational invariancc of the ground state. This 
provides the means to identify the ground-state phase 
boundaries. It had been rcalised^i^ that choosing interac- 
tion parameters depending on N oi L keeps the regions 
of interest centered, see figure [3] and figure 21 For the 
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FIG. 3: Results of exact numerical diagonalisation of the 
Bose-Hubbard Hamiltonian ([T]) with N = 5 bosons and — 1 
for various numbers of lattice sites L (order indicated by ar- 
rows holds for all panels) . The properties shown are indicators 
of qualitative changes in the ground state (cf.— i^i^): (top to 
bottom) the ground-state overlap with the non-interacting 
ground-state |x = 0), the incremental ground-state overlap 
(x|x + ^) (for A — 10~^), and the first excited energy relative 
to the ground-state energy. For explanations of the unitless 
axis in the middle graph see text. In this particular parametri- 
sation the transitional behaviour at Xci ~ 0.9, predicted by 
the non-linear Schrodinger equation approximation dicussed 
in Section Till Dl is visible. 




weak <~ attractive interaction T ~> strong 

FIG. 4: The same data as in figure [31 parametrised in terms 
of T. The dotted lines are guide to the eye, to mark three 
different regimes, the order of parameter L indicated by ar- 
rows holds for all three panels. In this parametrisation the L- 
dependence of the pre-transition coupling Td is apparent, as 
indicated by the thin dotted vertical lines. The pre-transition 
coupling Tc2 ~ 0.73, as indicated by the single thick dotted 
vertical line, is independent of L (cf.— ). At this coupling we 
see a minimum of the incremental ground-state overlaps and 
the onset of quasi-degeneracy of the ground and first energy 
levels. The numerical value is not in close agreement with the 
predicted value Tc2 = 2/3 of lllTEl This discrepancy may be 
explained by the fact that the particle number here is A'^ = 5, 
while the analysis in section UlI El assumes large particle num- 
ber. 



study of crossovers in the ground state the over-all en- 
ergy scale can be neglected, and we introduce parametri- 
sations mapping the whole region from the weak to the 
strong coupling limit into the finite interval [0, 1]. As a 
dimcnsionlcss coupling parameter of the model we de- 
fine (5 = 7/t, to study the ground-state properties of the 
model as 6 is varied. To help cope with the different scal- 
ing of the regions of interest as seen above, we introduce 
two further parametrisations in terms of dimcnsionlcss 
variables x, r S [0, 1]. These are defined by 

t = eAl-x), 7= If, (3) 

i = ^r(l-r), 7=^ (4) 

where e^, provide the energy scale. In terms of S we 
have 

NLS _ N5 

^~ 1 + NL6' ^ ^ 1 + N5' 

The non-interacting case is given by r = x = while 



pure interaction and no kinetic (hopping) contribution 
corresponds to t = x = 1- Other parametrisations, e.g. 
logarithmic dependence, are also used in the literature^^. 
Sec figure [3] and figure |4] displaying the same informa- 
tion, for a visualisation of the effect of the parametrisa- 
tion. Numerical exploration for small systems finds that 
the dips (local minima in the incremental ground-state 
overlap) arc quasi-stationary for scalings of 7 '-^ and 
of 7 ~ -Jj-, respectively. 



III. MEAN-FIELD THEORY 

Owing to the difficulties in treating the full quantum 
system with numerical and exact methods, approach- 
ing the system in the spirit of mean-field theory has 
been very popular. In particular with regards to in- 
vestigating non-linear phenomena like solitons, and de- 
scribing realistic experiments on BECs, these systems 
have been well studied in a wide range of contexts. The 
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continuum limit, known as the Gross-Pitaevskii equa- 
tionSS, the Lieb-Lingcr Bose gas or simply the non-linear 
Schrodinger equation (NLSE), have found wide inter- 
estii^. An extensive discussion of the mathematics of 
solution and further references for the discrete NLSE can 
be found in~. For the purpose of this paper we arc solely 
focused on pre-transitions in the ground state, though, 
and we will not consider these applications here. 

As the full discrete model is not integrable we will con- 
sider the three cases of the dimer (i = 2), trimer {L = 3) 
and the continuum limit {L s- oo). We will then com- 
pare these special cases with numerical solutions to the 
discrete mean- field equations for generic lattice size L. 
In the last part of this section we will present a semi- 
classical analysis, following a different approach^^. We 
will see that it recognises the second pre-transition not 
visible in the continuum limit, at least qualitatively, i.e. 
the critical interaction scales correctly with ^ jj, com- 
pared with in the continuum case. 



A. Generic L and A*' 

Consider the Hcisenbcrg equations of motion for the 
annihilation operators aj in ([1]) 



da, 



[a,,H] 



J 



I...L 



with restrictions to stationary solutions aj{t) = 
exp{—iEt)aj for some energy eigenvalue E. We make 
the usual mean-field approximation, here expressed by 
replacing the operators by complex numbers 



(operators) 



(complex numbers) 



The resulting -I- 1 coupled equations in the A^ -I- 1 vari- 
ables ai...aN, E for the finite lattice arc then given by 



E Qj = -t(aj_i + Oj+i) - 27|ajpaj , j = 

L 



(5) 



Note that we will discuss this procedure in Section fill Dl 
for the continuum model again. Either way, taking the 
mean-field approximation first and then going to the con- 
tinuum, or alternatively taking the continuum limit to 
the quantum Lieb-Liniger gas and afterwards replace op- 
erators by complex numbers, the result is the same con- 
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FIG. 5: Ground-state wave function overlap within the 
mean-field model of section IIII Dl versus attractive interac- 
tion strengh r. The numerical solution of the discrete non 
linear Schrodinger equation (solid black lines) exhibits the 
two dips locating the two transitional points, already antici- 
pated by figure[l]and figure[2] The continuum approximation 
via the exact soliton solution (section IIIID|I of the integrable 
non-linear Schrodinger equation (dashed red lines) has only 
one transitional point: this finite lattice effect is found in both 
the full quantum model and its semiclassical mean-field ap- 
proximation when comparing with each respective continuum 
limit. Note that the exact mean-field theory result does drop 
off again for strong interactions, which is not shown here. 
Confer figure [TT] for the exact solution in the case N = 2 
(Haldane-Choy). For strong interaction in this approxima- 
tion the ground state does not enter a region of small changes 
again, thus not specifying a second transitional point, nor 
does it relate to the location of the lattice model transitional 
point, see graph for L — 10. 



tinuum Gross-Pitaevskii equation. 

We show a numerical solution of these equations in 
figure [1] and figure O Clearly the limiting case of 
(de)localisation in the densities can be seen, as well as one 
distinct and one less sharp crossover within this mean- 
field theory. A clearer picture of the pre-transitions in the 
semi-classical analysis is given by the indicator property 
shown in figure El The two dips scale the same, namely 
~ jj and ~ respectively, found from exact diagonal- 
isation. As the discrete system ([S|) is non-integrable the 
general solutions are not known except for special cases. 
In the limit of weak interactions the ground state of the 
orginal system is given by the state with all particles in 
the zero momentum mode (an "ideal BEC"), correspond- 
ing to the constant solution = y^N/L, i ~ 1...N, 
with energy E = 2{t — 1 — t / L) (For this and the fol- 
lowing section we set er = 1). This delocalised wave 
function is a solution to the mean-field system for all 
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interaction strengths, but for stronger interaction the 
ground state becomes a locahsed solution. Higher energy 
solutions can be constructed by considering extensions 
Uj = yjN jL exp(27rj(/)/L) or sawtooth like amplitudes. 
Here we arc only interested in the lowest energy solution 
= 0. Although the particle number iV enters the sys- 
tem of equations as a parameter, this mean-field approxi- 
mation shows "sharp" pre-transitions between phases re- 
gardless of the value of N . In the following we will first 
discuss the special cases of the dimer (L = 2) and the 
trimer [L = 3), before considering the case of large L. 
This non-linear system has many solutions for a given 
parameter r. The numerical solutions shown were ob- 
tained by starting from a known limiting case and then 
iterating via small changes in r. As will be seen this pro- 
cedure may lead to spontaneous "hopping" to another 
solution once the current one ceases to exist. 



B. Dimer N = 2 

For the dimer system the mean-field equations ([5]) con- 
sist of three coupled, non-linear equations for the com- 
plex variables ai, ai and the energy E. Assuming real 
solutions is equivalent to finding the roots of a 4th order 
polynomial. Two solutions are real for the whole interac- 
tion range < r < 1, these lead to a\ = ±02, with "+" 
being the symmetric ground-state solution (the horizon- 
tal line in figure[B]). But for values of the coupling r > Tc2 
there opens up two new solutions with (the same) lower 
energy. These connect at Tc2 to the constant solution 

01 =02. For T — > 1 the solution localises, i.e. ai — > 1, 

02 ^ or oi ^ 0, 02 1, as shown in figure [S] The 
critical value Tc2 agrees with the semi-classical result of 
Sect. IHIEI and an alternative mean-field treatment given 



C. Trimer L = 3 

The trimer system is non-integrable and it has 
been studied previously in the context of chaotic be- 
haviour— Here we are only interested in soliton so- 
lutions for the ground state within the mean-field de- 
scription of the discrete non-linear Schrodinger equation. 
It is useful to introduce the notion of bright and dark 
solitons. A bright soliton has a localisation with pos- 
itive amplitude relative to the constant solution, while 
the dark solition has a negative amplitude, i.e. a bright 
soliton looks like a hill and a dark soliton looks like a 
valley. For the dimer case the twice degenerate bright 





Dimer 

1 . 1 . 1 . 1 . 1 . 1 . 1 . 1 . 1 


Trimer K 













0.2 0.4 0.6 0.8 1 



weak <— attractive interaction T — > strong 

FIG. 6: Occupancy of the modes for varying attractive inter- 
action r within the discrete non-linear Schrodinger equation 
approximation of section HIT AI Here the upper (lower) line 
shows either ai = az (a2 = as). The horizontal line de- 
notes the constant ("ideal BEG") solution, which exists for 
arbitrary L and r. In the dimer there exists a pre-transition 
coupling Tc2 = 2/3, beyond which a second, increasingly lo- 
calised solution exists. For the trimer case the pre-transition 
coupling in a two-height scenario is Tc2 ~ 0.663 (left dashed 
line). Note that the initially dark soliton (black lines) turns 
into a bright soliton before it merges with the (lowest lying) 
bright soliton (blue lines), cf. section Fill Gl for details. At 
Tc2 these real solutions cease to exist, but there is no smooth 
connection to the constant solution as in the dimer case, in- 
dicating the lattice effect, see text for further discussion, and 
cf. figure [2] 



soliton solution is at the same time a dark soliton, as the 
hill and valley cannot be distinguished for L = 2. For the 
trimer case L = 3 bright and dark solitons have different 
energies. 

Using a similar ansatz as for the dimer, i.e. 02 — 03 
and requiring ai to be real, reduces the problem to the 
analysis of a 4th order polynomial. A dark soliton (|a2| > 
|ai|) and a bright soliton (|a2| < |ai|) exist for r > Ta w 
0.663, with the bright soliton having the lowest energy. 
The energies of the two solitons are the same at t = r;, « 
0.692, which is not the point at which the bright and dark 
degenerate into the constant solution, i.e. |ai| ^ \a2\ at 
T ~ Ta- We remark further that at r = r;, the dark 
soliton becomes a second, higher energy, bright soliton. 
At this point the energies for this soliton and the constant 
solution are the same, as shown in figure [71 This ceasing 
of the real solution of the form ai > 02 = 03 is a hint to 
the qualitative difference between the dimer case L = 2, 
and the general case L > 2. We expect that in a small 
region the ground state is neither of the constant nor of 
the simple two-heights soliton form, but a more complex 
solution connecting these both. This intermediate region 
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FIG. 7: The energy for the lowest lying state of the two-height 
mean-field solution for dimer and trimer, corresponding to 
figure |6l Note that for L = 2 the soliton solution connects 
smoothly to the constant solution. In the two-height approx- 
imation for the trimer there is already a small region around 
T « 0.663 where the true ground state is not of the two-height 
form. Compare this to the large middle region visible in fig- 
ure [1] and figure [2l which differs from the conclusions in a 
recent study^. 



Schrodinger equation, where the latter reads 

Hnls^ f [*'t(x)*'(x)]dx 
Jo 

-c I '^\x)'^\x)'^{x)'^{x)Ax (6) 

^0 

with the periodic boundary condition ^'(0) = ^'(0- 
At this point we remark that the Hamiltonian ([5]) is 
integrable^^ - see section IIV CI for more details. One of 
the conserved operators is the total particle number 

N= ( ■^Hx)-^{x)Ax (7) 
Jo 

which is quantised and has eigenvalues which are non- 
ncgativc integers. The approximation of the Bosc- 
Hubbard Hamiltonian by the non-linear Schrodinger 
equation is 

Hbh ~ tA'HNLS-2W 

N ^ Af 



does not exist for the dimer. 



D. L ^ oo: Non-linear Schrodinger equation 
approximation 



In the limit as L — > oo we can approximate the Bose- 
Hubbard Hamiltonian ([T]) by a quantum field theory, with 
field operator '^{x) satisfying 

[^{x),^Hy)]^S{x~y). 



Setting aj ~ vA^E'Q -A), this consists of replacing 

— > J^^ dx under the assumption that A <C 1. 
This is to be understood as choosing L very large and 
N finite, distinct from the usual notion of the thermo- 
dynamic limit where N, L oo while keeping N/ L ~ 
constant. The implication of this approximation will 
be discussed later. These considerations lead to a map- 
ping of the Bose-Hubbard Hamiltonian to the non-linear 



where / = AL and c = 7/A<. Hereafter we set ? = 1 or 
equivalently A = L^^. 

The time evolution of the field operator "if can be de- 
termined in the usual way : 



dx'^ 



(8) 



Our next step is to treat ([5]) as a classical field equation 
(cfii). We introduce the rescaled field <I> = V N^^'^ and 
look for stationary solutions ^{x^t) = cxp(— ii?t)(f>(x) 
such that 



1 = / |$pda; 



^2$ 



dx"^ 



(9) 
(10) 



where U — 2cN. The ground-state symmetry breaking 
solution to equations (|9ll0p is knowniii^, and reads 



for U < 2tt^ 



$(a;) = 



lK{m) 
E{m) 



dn[2i4:(m)(a;-a;o)|m] for U > 2tt^ 
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Above, dn(u|TO) is a Jacobi elliptic function, E{m) and 
Kirn) denote the complete elliptic integrals of the first 
and second kind, and m is a function of t7— lii. Note 
that Xq G [0, 1] is the coordinate of the maximum of the 
wave function: the spontaneous symmetry breaking in 
the mean-field result is visible from the degeneracy of the 
solutions beyond the point of "collapse" of the constant 
solution into a soliton at the critical coupling Uc = 27r^. 
In terms of the dimcnsionless coupling parameter x of 
the Bose-Hubbard model this corresponds to 



5ci 



'nl 



or equivalently 



Tcl = 



Xcl = 



1 + 7r^ 



0.908 



(11) 



(12) 



The soliton solution connects continuously, but not 
smoothly, to the constant solution at Xc - in figure [5] 
this corresponds to dip in the dashed red line. A numer- 
ical solution for the finite size discrete NLSE is shown 
figure [T] and figure [21 here this corresponds to the sharp 
change at around r w 0.2. 



E. 



iV ■ 



oo: Semi-classical analysis 



minimise it subject to the particle number constraint 

L 



(14) 



The minimum occurs when 9j = O^j, which leads us to 
studying 



where 



It can be verified that for Nj = N/L\/j, Hi is glob- 
ally minimal and IHI2 is globally maximal. Thus for any 
7/i the solution Nj = N/LVj provides a fixed point of 
the system, which will be the unique global minimum 
when 7/i is sufficiently small. We look to determine the 
coupling at which this solution ceases to be the mini- 
mum. The results of the previous section indicate that 
when this happens a soliton solution will emerge. We can 
parametrise such a soliton solution as 



N 

Nj > — for j < z, 
ij 

N 

Nj < — for j > z 



(15) 
(16) 



In this section we present an alternative type of mean- 
field analysis, where we start by assuming that N is ar- 
bitrarily large and L is fixed. This is achieved by first 
canonically transforming to a number-phase representa- 
tion of the quantum variables. 

Let {Nj,dj}j=i,,,L, obey canonical relations [9i, 6j] = 



[N,, N,] = 0, [iV„ Ok] 
variables 

bj = exTp{i6j) 



5jk- We make a change of 



'Nj exp{—i6j) 



Using the fact exp{i6j)Nj = {Nj + l)exp{i6j) it can 
be verified that the canonical commutation relations 
amongst the boson operators bj, b'j are preserved. For 
large Nj we can approximate the Bose-Hubbard Hamil- 
tonian ^ by 



= -2tY^ yjNjNj+i cosiOj - Oj+i) 



(13) 



where 1 < z < {L — 1). Within this classical treatment, 
we can approximate the ground state for the full system 
by the two ground-state configurations for the sublattices 
j < z and j > z. For the full system at the pre-transition 
coupling Sci as predicted in llllDi we see that the systems 
on the sublattices are below the pre-transition coupling 
due to the L-dependence of 6ci ■ Hence the ground-state 
configuration across each sublattice is one where the Nj 
are constant on each sublattice. This leads us to look 
for soliton solutions within a two-height approximation, 
valid close to a point of broken symmetry: 



N{l + a) 

N{1 - a) 
2{L~ z) 



for j < z 
for j > z 



where — 1 < a < 1 is continuous, such that holds. 



We now treat H as a classical Hamiltonian and look to 



In terms of the above parametrisation the Hamiltonian 
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IS 



= -tN 



(z-l)(l + a) (L-z-l)(l-a) 
z L — z 



1-^2 



(17) 



y z{L-z) I 

-fN^ fil_+af_ (1 -a)2 
z L — z 



The next step is to minimise this expression for H with 
respect to the variables z and a, this involves mostly 
standard calculus techniques. We find that the small- 
est coupling for which symmetry-breaking of the ground 
state occurs is 



5c2 = 



2L 



Xc2 



1 + 2L' 



Tc2 



(18) 



In deriving the value Sc2, we justified the use of the two- 
height soliton approximation on the basis that Sd is a 
decreasing function of L. The fact that Sc2 was ultimately 
found to be independent of L does not invalidate the two- 
height approximation within the context of the above 
analysis. The numerical results of figure 2] illustrate that 
in a generic finite system the gap between the ground 
state and first excited-state energy levels is smaller at 
Sc2 (or equivalently Tc2) than at Sci (or equivalently Td). 
Thus the notion of quasi-degeneracy of the energy levels 
is more appropriate at dc2 than at 6ci- 

Let us consider what happens when we now take the 
thermodynamic limit N, L oo: 



Sci = 0, Xci 



l + 7r2' 



Tcl 0, 



r.2 - 3. 



(19) 
(20) 



Sc2 = 0, Xc2 = 1, 



The fact that the two sets of values (|T9l) and ([20|) do 
not agree is an indication that the limits N 00 and 
L — !■ 00 do not commute, meaning that the usual con- 
cept of the thermodynamic limit is not well-defined for 
this model. Equations (fT9|) and ([20]) again show that 
two of the regions of interest will vanish when using the 
standard 6 variables. When using the parametrisation 
in the variables x and r only one region disappears and 
one pre-transition point stays finite, i.e. away from (no 
interaction limit) and 1 (no hopping limit). 

Another curious point to observe is that while both Xc2 
and Tcl are L-dependent, they are in fact independent of 
N. This gives faith that the general qualitative ground- 
state features of the finite system will be tractable from 



q-Oscillator 



Asymmetric 
discretization 



Symmetric 
discretization 



integrable extentions I for general N,L 



t 



Bose Hubbard model 

non-integrable for general N,L 



Haldane-Choy N=2 



FIG. 8; Integrable relatives of the Bose-Hubbard model: the 
6 small boxes have Bethe Ansatz solutions, while the general 
Bose-Hubbard models is non-integrable. 



analyses of systems with relatively small particle num- 
bers. 



IV. LIMITING INTEGRABLE MODELS 

While repulsive boson systems, as well as repulsive 
and attractive fermion systems are well studied in the 
context of solvable systems the attractive boson gas re- 
ceived comparatively little attentionii. Still the seminal 
Bethe Ansatz solution for one-dimensional contact inter- 
action bosonsi^ in the continuum describes the repulsive 
as well as the attractive regime, in which the solutions to 
the Bethe Ansatz equations become of different character 
Initially it was believed that also the Bose-Hubbard 
model has a Bethe Ansatz solution^^, but it soon turned 
out that this model is non-integrable. Nevertheless there 
are several integrable limits and extensions, see figure [H 
Out of these we will examine the three limits shown inside 
the general Bose-Hubbard box for the case of attractive 
interactions. Integrable lattice distortions of bosons on a 
one-dimensional lattice, for instance the three boxes on 
top of figure[Sl have been studied mainly for the repulsive 
(,g^gg25j30^^i46_ fpj^g attractive parameter region is tech- 
nically harder than the repulsive case: for instance the 
attractive Bethe Ansatz roots in the Licb-Liniger model 
lack several of the properties which allowed analysis for 
repulsive interaction: e.g. string solutions which keep 
their string form and saturation of the root distribution 
in the thcrmodynamical limit. The problem of collapse of 
the system already at infinitely weak attractive interac- 
tion when taking the thermodynamic limit is less prob- 
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lematioi^ and has been addressed by the iV-dependent 
reparametrisationi^. In the next secion IIV Al we will 
briefly discuss the Bethe ansatz solution of the L = 2 
dimer case and point out that it has only one crossover. 
Then we move on to discuss the finite lattice size case 
L > 2 via the two-particle solution of the Haldanc-Choy 
ansatz: here we see now two pre-transitions, i.e. dips in 
the indicator property ground-state overlap, which have 
the correct scaling behaviour found numerically and via 
semi-classical analysis earlier in this paper. 

We will establish via the N = 2 Haldane-Choy solution 
and the results from section IIIII for the relation between 
the (discrete) NLSE and the (continuum) GPE that the 
integrable Licb-Linigcr continuum gas is a good proxy 
for the discrete lattice model for the study of the NL- 
dcpcndcnt pre-transition. This motivates our discussion 
of the Bethe ansatz root distribution for the attractive 
ground state of the Lieb-Liniger model and relates these 
results to the Bose-Hubbard model in the last section. 



Bose-Hubbard dimer 



erator 

H = -2t(7Vw+ (1-1*2)4^ 
\ au 

d^ d 
-7(2^2 + 2(1 - N)u— +N'^ -N) 
au^ du 

d2 d 
= -27^2— + (27(iV - l)u + 2i(u2 - 1)) 

dw^ du 

+ (iV7 - N'^-i - 2tNu) . (21) 

Now we look for solutions of the eigenvalue equation 

HQ = EQ (22) 

where Q is a polynomial function of order N which we 
express in terms of its roots {vj}: 

N 

Evaluating (|22p at u ~ Vk for each k leads to the set of 
Bethe ansatz equations 



For the dimer case, L = 2, the Hamiltonian ([T]) reduces 

to 

H = -2t(^a\a2 + a^aij — 7 (^a\a\aiai + a\a\a2a2 

This Hamiltonian can be expressedi^ in terms the su{2) 
algebra with generators {5^, S^} and relations 

[S\S^] ^ ±5*^, [S+,S-] = 2S'. 

Using the Jordan-Schwinger representation 

= a\a2, S~ = a\ai, S^' ~ ~ ^^2^2) 

this leads to 

H = -2t {S+ + S-)-j (^2{S' f + ^N^ - . 

The same (iV-|- l)-dimcnsional representation of su{2) is 
given by the mapping to differential operators 

d iV 2 d d 

S ^ u- — , b —Nu — u — , b ~ — 



N 



Vj - Vk 



By considering the terms of order N in ((22|) the energy 
eigenvalues are found to be 



N 



E = jN(l- N) + 2tY^ 
i=i 



(24) 



We can transform the differential operator (|2T|) into a 
Schrodingcr operator—. Setting 

* cxp cosh(y27a::) - ^ ^j^^ 



N 



du 2 



du 



du 



where the potential V{x) is 



acting on the space of polynomials with basis 
{1, u, w^, u^}. We can then equivalcntly represent the 
dimer Hamiltonian H as the second-order differential op- 



V{x) = ( ^ sinh2(y2^x) - 2t{N + 1) cosh(v/2^a:) j , 
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O 0.4 

V 



N=400 



T 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 

weak <— attractive interaction x - 



0.8 0.9 
> strong 



FIG. 9; The Bose-Hubbard Dimer L — 2 has a single 
pre-transition point. Shown are (top to bottom), ground- 
state wave function overlap with the non-interacting reference 
state, the incremental ground-state wavefunction overlap and 
the first excitation energy relative to the ground-state en- 
ergy, each obtained by exact numerical diagonalisation. The 
dashed line indicates the theoretical value of Tcrit =2/3 which 
is given by mean-field theory. For the wavefunction overlaps 
in the middle graphic the mean-field result for the overlap 
are shown by the dotted line, which is barely distinguishable 
from the numerical result. The value obtained from the exact 
Bethe ansatz solution is tba ~ 2/(3(1 + N~^)), giving the 
quantum correction to the mean-field result. 



then 



dominate, while after a smaU crossover region for large 
T this non-interacting BEC state has a very low relative 
weight in the ground state. In the complementary plot, 
against the iV-body Schodinger cat-state |iV, 0) -I- |0,iV) 
as reference state, the overlap would be almost constant, 
close to 1 in the strong interacting regime on the right, 
while it would be ^ 1 and quasi-constant in the weak 
interacting region: the relative weight of a localisation 
of N particles is low. Similar for the bottom picture in 
figure [9] the ground-state energy for the very weakly 
interacting regime r « is non-degenerate, the first ex- 
citation is separated by the energy required to transfer a 
single boson from the zero momentum mode to the first 
momentum. In the strong interacting limit for large r the 
ground state is quasi degenerate : the (anti-)symmetric 
cat states |iV, 0) ± |0, A'') have the same energy. 

The mean-field calculation for the dimer, see figure [6l 
shows the (square root of) the relative occupancy of 
the two sites. For r below the critical interaction both 
sites are equally occupied, the totally delocahsed con- 
stant solution, with all the particle in the lowest mo- 
mentum mode bo- At Tc = | the symmetry breaks and 
one site has higher occpuation than the others. Due to 
the quantum-mechanical superposition in eigenstates this 
can only be seen in the mean-field theory. The second 
momentum mode has a finite and increasing occupation 
beyond the critical interaction, though, and it reaches 
?ifc=o = nk=i = i for r ^ 1. This is the complete de- 
localisation in momentum space and corresponds to the 
complete localisation in real space density observed in 
the soliton solution. 



with E given by (j24p whenever the {vj} are solutions 



of the Bethe Ansatz equations (p3l) . It is easily checked 
that the potential has a single minimum when j /t = S < 
2/{N + 1) and two minima when 6 > 2/{N +1). The 
critical value 6 = 2/(A^+ 1) agrees, to leading order in N, 
with the mean-field theory result for Sc2 of section IIIIEl 
cf. equations (fT8|) . 

From the analysis of the Bethe Ansatz solution it can 
be seen the limiting case of only two lattice sites has a 
single transitional point, visualised in figure [S] This cor- 
responds to the case where the two minima in figure [3] 
coincide and the middle region is no longer visible. In fig- 
ure [5] two physically very different regimes can be seen: 
the ground-state overlap measures the relative weight the 
occupation of the zero momentum mode by all particles 
has (relative to the non-interaction BEC state with 100% 
condensation at t = 0). For small r these contributions 



B. Haldane-Choy Bethe Ansatz for = 2 

The Bose-Hubbard model ([T]) has a Bethe Ansatz so- 
lution in the spirit of the fermionic Hubbard model, but 
it is only solvable for a maximum site occupation of two 
particles^S. For iV = 2 the exact eigenstates are 



|BA) 



L 



^i{kn+qm) , sin fc-sin q-if ^i(Qn+km) n < m 
sinfc— sing+i7 ' — 

, n > m 



with the Bethe Ansatz equations 



^■^^ ^ sinfc-sing-i7 ^ ^^^^ ^ sing-sinfc-i7 ^^^^ 



sinfc — sin(7+i7 



sing— sin fc + i7 
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to i» 



L=10,50,100,1000,10000 (Bcthe Ansatz) 
L-10,50,100 (exact numerics) 



0.2 0.4 0.6 0.8 1 

weak <— attractive interaction T — > strong 

FIG. 10: Similar to figure|4]for the solvable case of two bosons 
in a lattice of differing size L (order indicated by arrows holds 
for all three panels) Already for the minimal particle num- 
ber N = 2 the two pre-transition points can be clearly seen, 
compared to the one pre-transition point for the dimer, see 
figure |9] The exact Bethe Ansatz solution ((28} (solid lines), 
available for arbitrary L, is compared with exact numerics 
(dots) for small sizes. Shown are (top to bottom): the ground- 
state overlap with the non-interacting ground state |t = 0), the 
incremental overlap {t\t + A) for A — 10~^, and the first ex- 
cited energy relative to the ground state L^{Ei — Eq). 



Here — a|^aj|0) is not normalised, i.e. («, j) = 1 

respectively — 2 for i j respectively i — j ■ The energy 
eigenvalues for these are E = — 2(cosA: -I- cosg), moti- 
vating the name "quasi-momenta" for the Bethe Ansatz 
roots k and q. The Bethe Ansatz roots for the ground 
state are symmetric, k = —q, and imaginary for attrac- 
tive interactions 7 > 0. For this setting we can define 
k = —q = iK , with K > determined by the single 
Bethe Ansatz equation, here in inverse form 



j{K) = 2sinhA' tanh 



KL 



(26) 



For use in the next section we also note that for 7 > 
^ cos ^ the two real roots of the BAE solution for the 
first excitation Ei merge to a complex 2-string of the form 
k,q = ^ ± iK with K > 0, and the inverse function given 
by ■'y{K) = 2 cos ^ coth sinh/\. The real roots of the 
first excitation for 7 < ^ cos ^ are given hy k = K, q — 
^ — K with < K < j^. The inverse function relating 
the parameter K to the interaction strength is 'y{K) ~ 
2 cos ^ tan sin ^^j^^ ■ We use these expressions for 



the analysis of the indicator properties like L^{Ei — Eq), 
in figure 1101 as well as for comparison with the Lieb- 
Linigcr continuum model in the next section. 



The (not normalised) ground-state wave function can 
be written as, cf. ((32|) : 



\K) Jeif(#-|n-m|) +e-/^(*-|n-™|)j ^27) 



E 



2 coshii'(|n — iti\ — —) 



\n, m) 



resulting in the closed form expression for the (not nor- 
malised) overlap in the Haldanc-Choy model 

(A' + A|A'-A)= (28) 
AL (coth K coth KL + coth A coth A L) 

together with the normalisation 

{K\K) = Lc-'^^ (coth A' - 1) 
X (2Lc^(^+2) - 2Le^^ + e^^'^+i) + e^^^ - e^^ - l) 

this results in the normalised overlap expression 
{K + A\K~A) 



+ A\K + A){K - A\K - A) 

In the above equations K ± A denote the imaginary 
parts of the single Bethe Ansatz root associated with 
the two different interaction strengths ti i—t K + A and 
T2 1-^ K — A, i.e. solutions to ([26|l . This expression 
depends on the interaction strength 7 only through the 
Bethe Ansatz root K, allowing closed form solution in 
parametrised form^. The Bethe Ansatz solution for only 
two particles is not truly a many-particle solution - the 
N = 2 Bose-Hubbard model can be treated exactly con- 
ventionally in center-of-mass coordinates^!^. In that 
case the physical meaning of the Bethe Ansatz quasi- 
momenta is lost, though. The solution presented here is 
visualised in figure 1101 Wc remark that within this ap- 
proach the exact momentum distribution of two bosons in 
the one-dimensional lattice can be calculated explicitly, 
clarifying the connection between the (here two) Bethe 
Ansatz quasi-momenta and the physical momenta, which 
is of interest for example in the intcgrable boson-fermion 
mixture^. 
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C. Lieb-Liniger approximation 

The continuum model in ^ is the integrablc Lieb- 
Linigcr gasi^. For the repulsive regime it is arguably 
one of the best studied integrable model a-^^'-^^i^^'^''i^^'^^ , 
while the attractive regime is less popular—, due to diffi- 
culties in taking the thermodynamic limit. When taking 
the limit L ^ oo in the Bose-Hubbard model the Lieb- 
Liniger model can be used as an integrable approximation 
for the weak coupling limit. Information is lost in taking 
this limit, i.e. in going from the three independent pa- 
rameters N, L,7 of the Bose-Hubbard model to the two 
parameter continuum model. Thus we expect the range 
of validity to be restricted, but in turn the property of 
integrability is gained. There are two different ways of 
looking at this integrable model as a limit of the Bose- 
Hubbard model, we consider the analysis for finite N and 
large L — > oo in section HV C II For repulsive interactions 
the thermodynamic limit can very succesfuUy be treated 
(see2^ for references) for constant density n — when 
N, L oo. In the first case the two independent parame- 
ters are A'' and an interaction strength. In the second case 
we keep the particle density/filling factor n and an inter- 
action strength as free parameters. This physical notion 
of density (of Bethe ansatz roots) is not extendable to 
the attractive case, the bosons tend to cluster up instead 
of saturating. We discuss this further in section [TV C 21 



1. Analysis for finite N and large L — > cxd 

In this section we analyze the special case iV = 2 as ex- 
ample for finite TV and (very much) larger L ^ N . The 
exact Haldane-Choy solution discussed in section IIV Bl 
is the yardstick to explore the impact of the continuum 
approximation in the full quantum model. The energy 
eigenvalues in the general Lieb-Liniger model correspond- 
ing to an approximated lattice model are given by 



N 

L^E = ^ fcf + const. 

i=l 



(29) 



with the N complex parameters ki determined by solu- 
tions to the Bethe ansatz equations 



N 



■j-j- kj kj ir-^ 

j^i ^ + ^'''iv 



i = l...N (30) 



In particular we see that the continuum model gets 
mapped onto the weak coupling limit of the lattice model. 



as ^ — > for N = 2 and L — > oo. To check how well 
the Lieb-Liniger model approximates the Bose-Hubbard 
model for large L we calculate analytically the ground- 
state overlap for the case N = 2 and compare with the 
Haldane-Choy expression The root behavior for 

ground state and first excitation (see also appendix ofi^) 
is similar to the lattice case: the two ground-state roots 
form again a purely imaginary complex pair fci 2 = iiX, 
where the inverse function is given by 



7 



K K 
2 — tanh — . 
L 2 



(31) 



The first excitation roots form a complex pair past the in- 
teraction strength 7 > 7c = of the form fci.2 = TriiX, 
with inverse function 7 = 2-^ coth . For weak inter- 
action 7 < 7c the first excitation has two real roots 
at k2 — 2-K-ki — 2-K-K^ with inverse function 7 = 
|;(7r-/C)tanf . 

The (not normalised) ground-state wave function for 
finite interaction and = 2 is given by, cf. 



=2coshA'(|a;-y| - i) 



(32) 



The normalised ground-state wavefunction overlap for 
two different interaction strengths, with corresponding 
imaginary part of roots ± A, is given by. cf. 



+ A|X- A) 
/sinh(/<) sinh(A)\ 



(33) 



i^2- A2 



y (A' - A smh[K - A)){K + A + sinh(A' + A)) 

The comparison for iV = 2 of the rescaled Lieb-Liniger 
gas with the Bose-Hubbard system is shown in figure [TT] 
From the exact diagonalisation of small systems, see fig- 
ure [3] and figure \W\ it is apparent that for increasing 
lattice size L and fixed particle number A'' a growing re- 
gion extends from the strong interaction limit r = 1 to 
smaller r. The shown physical properties in this region 
are independent of L, nevertheless the L — > 00 Lieb- 
Liniger model is not a valid approximation for that re- 
gion: from the Bethe Ansatz equations and ([50)1 it 
can be seen that the interaction strength gets rescaled 
by ^ , effectively mapping the Lieb-Liniger model 
onto the infinitely weak interacting Bose-Hubbard model 
by quasi- linearising the roots. The continuum limit does 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



weak <— attractive interaction x — > strong 

FIG. 11: Ground-state wave function overlap versus attrac- 
tive interaction strength r for N = 2 bosons. The exact 
solution of the full quantum model (black solid lines) on the 
finite lattice (Haldane-Ghoy) exhibits two minima, indicat- 
ing two transitional points. The continuum approximation 
via the Lieb-Liniger model (red crosses) discussed in sec- 
tion |IVCT] displays only one minimum, indicating it has only 
one transitional point, see text and cf. figure [5] for the large 
N mean-field result. Note that the agreement is best in the 
L-dependent weak regime (left side) and not in the the L- 
independent strong interaction regime (right side). 




root value K 



FIG. 12: Quasi-distribution of the Bethe ansatz roots n{K) 
for numerical solution of (|35|l with interaction F below the 
pre-transition value Fc = tt^. The roots follow the semi-circle 
law p7|) . found analytically in the weak coupling limit F — > 0. 
The numerical results for A'' = 51 are in agreement with the 
notion that the semi-circle law holds asymptotically for F < 
Fc for sufficiently large A'', while for F > Fc the distribution 
is uniform. The dashed blue lines show the start of the finite 
size crossover from the semi-circle towards the box shape (|38[). 
The inset is an enlargement of the center region where the 
inner-lying roots approach the uniform density first. 



not capture the physics in the strong interaction region, 
in particular it does not see the second pre-transition 
point Tc2 connected to the finite lattice efi'ect. 

Our results for the semi-classical and the exact solution 
in the two-particle sector for the quantum model indicate 
that this limit is useful for the first crossover, though. 



2. Analysis of Lieb-Liniger equations for large N 

The complicated form of the wave function within the 
coordinate Bethe Ansatz makes a straight forward ex- 
tension of the ground-state wave function overlap calcu- 
lation similar to equation p3|) for the two-particle case 
impossible. There exist determinant formulations via the 
Algebraic Bethe Ansatz^. These have been studied for 
the repulsive case only, though. 

In the general Lieb-Liniger equations the first excited 
energy is relatively complicated to treat, as the root pat- 
tern is not as simple as in the N — 2 case. It can be 
shown that the roots never merge into the true A^-string 
for N > 2 for total momentum one, as this would vio- 
late hermiticity^. The ground-state root configuration is 
more accessible, as it is generally believed to be an ideal 
A^-string: the roots are purely imaginary and distributed 
symmetrically around the origin. The limit of strong in- 



teraction or very large box size L has been previously 
studied: the roots are then asymptotically linear in the 
interaction^ and evenly spaced. In this limit the (not 
normalised) wave function is of the McGuire fornJ^i^ 



vI/(x)=exp(-M^ 



(34) 



which is also relevant to (infinite length) optical wave- 
guides^. 

To analyse the attractive ground state of an abstract 
system of equations of Lieb-Liniger type we introduce the 



real variables Kj via kj 



iKj . It is ad- hoc assumed that 



for finite interaction and finite N there exists a unique 
real solution to 



e^- = 



N 



T_ 

N 



T_ ' 

N 



= I...N . 



(35) 



Here T = cN is the rescaled interaction. Note the sim- 
ilarity of these ground state equations to systems with 
hard wall boundary conditions^i^^ due to the externally 
imposed symmetry of the roots. The formulation of the 
problem in terms of the variables {Ki\ and T allows the 
definition of a sensible distribution or quasi-dcnsity of 
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Bethe Ansatz roots for very large but finite particle num- 
bers N. Here it is useful to define a quasi-density of the 
roots Ki, which for example can be done via^ 



n{x) ^ { ^' 



-Ki 



1^1 ^ -^max 



(36) 



In the weak coupling limit the root distribution of the real 
solution to ((35|) follows a semi-circle law derived from the 
relation to the Hermitc polynomials^LM 



1 



~K\ \K\ < /w = 2VV (37) 



In the strong interaction limit the application of the 
string hypothesis leads to a uniform, box shaped den- 
sity. When constructing a string solution to for 
fixed N and increasing T oo the difference between 
closest roots is asymptotically i^i+i — Ki = j^. Sum- 
ming up over the symmetric root distribution it follows 
that -RTmax = i'n^ ~^ \^ this agrees with numerical 
exploration for small particle numbers iV < 50. 



n{K) 



1 



2Kr, 



\K\ < K,, 



(38) 



The above expansions hold in the limits of F ^ (weak) 
and r ^ oo (strong), respectively, while N (large) is held 
constant. Note that i^max is in both cases independent 
of the particle number N. This would allow, at least 
in principle, to explore the interesting limit N ^ oo for 
a fixed and finite interaction strength < F < oo. It 
is technically hard to relate the Bethe Ansatz roots in 
Lieb-Liniger type models to physical properties within 
the exact approach. Nevertheless, it is expected that the 
quasi-distribution of the Bethe Ansatz roots in (pS)) 
for large N (resp. N — > oo) will show qualitatively dif- 
ferent behaviour in the two regions F < Fc and F > Fc, 
see figure [T^ for numerical results for 51. 

For weak interaction F <C Fc the numerical solution 
{Kj} is distributed approximately as a semi-circle ((37)l . 
while for F ^ Fc the quasi-density approaches a uniform 
box shape ([55)1 . Numerical results for small system sizes 
suggest an agreement with the expected value Fc = tt^ 
separating the two regions, where Fc is the location of the 
single minima in the ground state wave function overlap 
in the continuum model as discussed in the earlier sec- 
tions of this paper. 

Numerical solutions to finite N Lieb-Liniger equations 
are usually found by starting with an initial guess of the 
roots in a known region, e.g. the weak coupling limit. 
Then the interaction is increased in small steps F —f F+A 



where N remains necessarily constant. This so called 
root tracking works well if the root set {fci}|r-i-A is sim- 
ilar to the previous step {fci}|r - for a close initial guess 
most non-linear solver have good convergence. From the 
above it can be seen that this method is unsuitable for the 
study of large N behavior - the root pattern is expected 
to change strongly when crossing over the pre-transition 
at Fc = TT^ , which is in accordance with findings of Sak- 
man et a\^. In a diagram A^ vs. F the above method 
corresponds to moving along horizontal lines, where in 
the left part the root distribution is asymptotically of 
semi-circle shape, while on the right side it has the uni- 
form box shape. Using that the quasi-density for 
finite particles is in one-to-one correspondence with the 
root set {Ki} the system can be solved on vertical 
lines, i.e. for fixed interaction F and increasing A^. In 
that way the solution is stable, i.e. it does not change 
significantly for increasing A^ as the pre-transition point 
is not crossed. 

The preliminary numerical results obtained by root 
tracking agree with the behavior described above. Nev- 
ertheless, a rigorous analysis of ([55]) is necessary to deter- 
mine if a quantum phase transition occurs. In particular 
the critical value Fc — tt^ has not been found from the 
Bethe Ansatz equations. This result will be relevant for 
the description of the first pre-transition of the initial 
finite size Bose-Hubbard model, when transforming the 
considered abstract system back to the physical problem. 



V. CONCLUSIONS 

In this paper we have argued that there are signs of 
transitional behavior in the ground state of the attrac- 
tive one-dimensional Bose-Hubbard model. A discussion 
using conventional Quantum Phase transitions - defined 
in the thermodynamic limit of many particles A^ on many 
sites L - is unsuitable as the standard limit for attractive 
bosons is subject to instant collapse. Instead we have 
used the notion of pre-transitions, characterised by a sud- 
den change in the ground-state properties when crossing 
a threshold interaction strength in a system of large but 
finite size A^, L. 

Such pre-transitions are visible in indicator properties 
as for instance the energy gap between ground state and 
1st excitation indicating onset of degeneracy, and local 
minima in the incremental overlap (r -|- A|t), where |t) 
is the ground state for attractive interaction strength r. 

We have used mean-field like approximations and 
intcgrable limits of the model to examine regions in- 
accessible to exact diagonalisation, and compared with 
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exact numerics where applicable. The transitional region 
depends on both lattice size L and number of bosons 
iV in a non-trivial way. For specific parametrisations 
of the coupling strength between the kinetic and the 
interaction contributions in the Hamiltonian one of 
the crossover points is quasi stationary while the other 
wanders. In particular we have shown that in the limit 
of very small and very large lattice size L, the complex 
transitional regime reduces to only two regimes with 
one single crossover point, in agreement with earlier 
studies on these models. The ground state is predicted 
to change strongly in a small region around critical 
attractive interactions. In experiments with controlled 
change of attractive interaction this should have clearly 
visible effects in properties like correlation functions and 
momentum distribution. If ultracold quantum gases with 
large but finite particle number TV and lattice size L, 
enter the strong attractive interaction region the validity 
of the physical description by the simple Bose-Hubbard 
model needs to be carefully investigated, though. In 



addition it will be interesting to see how this transitional 
behavior manifests in theories of more complex attractive 
boson systems, as we believe this is a generic feature of 
attractive bosonic systems rather than a speciality of this 
particular model. The generalisations already studied in 
the repulsive regime like long-range hopping, long-range 
interactions and extensions of lattice geometry to ladders 
and square lattices are an obvious starting point for 
further exploration. For these systems there are cur- 
rently few methods available using integrable techniques. 
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